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Abstract 

Many non-Hermitian but PT-symmetric theories are known to have a real 
positive spectrum. Since the action is complex for there theories, Monte 
Carlo methods do not apply. In this paper the first field-theoretic method 
for numerical simulations of "PT-symmetric Hamiltonians is presented. The 
method is the complex Langevin equation, which has been used previously to 
study complex Hamiltonians in statistical physics and in Minkowski space. We 
compute the equal-time one-point and two-point Green's functions in zero and 
one dimension, where comparisons to known results can be made. The method 
should also be applicable in four-dimensional space-time. Our approach may 
also give insight into how to formulate a probabilistic interpretation of VT- 
symmetric theories. 
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I. INTRODUCTION 



Traditionally, only theories with Hermitian Hamiltonians are studied in quantum me- 
chanics and quantum field theory. This is because Hermiticity guarantee real eigenvalues 
and therefore unitary time translation and conservation of probability. It has recently been 
observed that quantum mechanical theories whose Hamiltonians are not Hermitian but are 
symmetric under a transformation known as PT symmetry have positive definite spectra 
||T|-rB|. A major criticism of these theories is that a consistent probabilistic interpretation 



has not been formulated. This paper suggests that there is a real, Fokker-Planck probabil- 
ity underlying these theories and presents a numerical method for calculating the /c-point 
Green's functions, Gk-, of these theories. 

A PT-symmetric Lagrangian that has been studied in the past is defined by the Euclidean 
Lagrangian 

C, = \m' + \m^^-j^m^. (1.1) 

A recent paper used Schwinger-Dyson techniques to study this self-interacting scalar quan- 
tum field theory [|1^. Green's functions, G^, calculated by this method agreed extremely well 



with known results. It was argued that these theories possess a positive definite spectrum 
and a nonvanishing value of Gi = (O|0|O) for all N > 2. 

Under VT symmetry V sends — —0 and T sends t ^ —t and i ^ —i, where t is 
time. (Note that in one dimension (i.e. quantum mechanics) represents the position of the 
particle and V corresponds to refiection in space.) Thus, Ce is manifestly VT symmetric. 
It is believed that the reality and positivity of the spectra are a direct consequence of this 
VT symmetry. The positivity of the spectra for all is an extremely surprising result; it is 
not at all obvious, for example, that the Lagrangian Ce = {d(j)Y /2 — gcf)^ / A corresponding to 
A^ = 4 and m = has a positive spectrum. To understand this and other results, we must 
properly define the contours of integration for path integrals and, in dimension D > 0, the 
boundary conditions in the corresponding Schrodinger equation. 
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Contours of integration and boundary conditions have been extensively studied in zero 
and one dimension [0,0]. In zero dimensions one must choose the contour of integration for 
the path integral such that it lies in a region where exp(— 5'(0)) is damped as |0| — > ±00 
and the path integral converges. For massless versions of Eq. ( |1.1[ ), these regions are wedges 
and are chosen to be analytic continuations of the wedges for the harmonic oscillator, which 
are centered about the negative and positive real axes and have angular opening ^. For 
arbitrary N > 2 the anti-Stokes' lines at the centers of the left and right wedges lie below 
the real axis at the angles 

'(iV-2)7r\ 



6'lcft = -TT + 



2N 



^^right - - 2iV '* ) ■ '-'^■^^ 

The opening angle of these wedges is 

Similarly, for one-dimensional versions of Eq. (|1 . 1| ) with m = 0, the Schrodinger differ- 
ential equation is 

In Ref. it was shown how to continue analytically in the parameter N away from the 

harmonic oscillator value N = 2. This analytic continuation defines the boundary conditions 

in the complex-0 plane. The regions in the cut complex-0 plane in which ip{(f)) vanishes 

exponentially as |0| —>■ 00 are once again wedges. These wedges also define the regions 

in which exp{—S {</))) is exponentially damped and the corresponding path integrals are 

convergent in one dimension. Once again, the wedges for > 2 were chosen to be analytic 

continuations of the wedges for the harmonic oscillator. For arbitrary N > 2 the anti-Stokes' 

lines at the centers of the left and right wedges lie below the real axis at the angles 

fN-2\7c 
ViV + 2y 2' 

(1.4) 

The opening angle of these wedges is 
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N + 2J 2 



Consequently, expectation values for PT-symmetric theories can be understood as path 
integrals that have been analytically continued in A^. This analytic continuation deforms the 
contour from the real axis for the harmonic oscillator, = 2, to contours in the complex- 
(f) whose endpoints lie in wedges where exp(— 5'((/))) is damped as \(j)\ oo and the path 
integral converges. Defining the complex variable (pc to follow any contour whose endpoints 
lie in the appropriate wedges, PT-symmetric expectation values of operators A = A{(j)) are 
given by: 

(0|^|0) ^ ID<l>c A{<Pc) e-^('^-) 
(0|0) JD(f)c e-^i^c) ' ^ 

where S{(f>c) = J d^XCE[(pc{^)] = I d^XH[(f)c{X)] is the Euclidean space action. The 
most common choice of contour is the one along which exp(— S'(0)) is purely damped. This 
contour is defined as (pc = exp(i^^/^), — oo < r < and (pc = exp(i^^/j), < r < oo, where 
Ol and Or are defined by Eq. or Eq.(|r|). 

As mentioned, another remarkable property of the Lagrangian in Eq. ( |1 . 1| ) is that for 
all A^ > 2 the expectation value Gi = (O|0|O) of the position operator in the ground 
state is nonzero. This surprising result shows that the theory is not parity symmetric even 
when A^ is even. The violation of parity symmetry is a consequence of the manner in which 
the boundary conditions and path integral contours are defined. The boundary conditions 
require that as \(f)\ —>■ oo the anti-Stokes' lines are in the lower half of the complex-0 plane 
for N > 2. While Eq. ( |1.1|) appears to be invariant under a parity transformation, the anti- 
Stokes' lines are sent to the upper half of the complex-0 plane; this corresponds to a different 
set of boundary conditions. Thus, the theory is not parity symmetric except for the special 
case of A^ = 2 where the anti-Stokes' lines are the real axis. Notice that if we now perform 
a time reversal transformation, complex conjugation sends the anti-Stokes' lines back down 
into the complex-0 plane, and both Eq. ( p. . 1|) and the boundary conditions are identical to 
the original formulation of the theory. That is, these theories are PT-symmetric, but they 
are not symmetric under V or T separately. As a result, Gi is pure imaginary and G2 is 
real. 
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The results for PT-symmetric but non-Hermitian theories suggest that these completely 
new theories may describe physical processes. Previous studies have obtained extensive re- 
sults in zero and one dimension (i.e. quantum mechanics) but have been unable to perform 
calculations in higher dimensions. Here, we use the complex Langevin method to calculate 
the equal-time one-point and two-point Green's functions for massless versions of Eq. ( |1 . 1[ ) 
in zero and one dimension with = 3 and = 4. The results are in good agreements 



with those computed by numerical integration [|I| and by variational methods |114| , p!7[| . Ref- 
erence studies complex Hamiltonians using the Langevin approach in higher dimensions 
and also obtains accurate results. This suggests that the complex Langevin method is the 
robust numerical method needed for studying PT-symmetric theories in higher dimensions. 
Consequently, we are currently performing field-theoretic, numerical simulations for VT- 
symmetric Hamiltonians in two space-time dimensions. The current paper also reveals an 
underlying real Fokker-Planck probability for these theories. We believe our results represent 
a significant step towards a physical understanding of PT-symmetric theories. 

This paper is organized as follows: In Sec. y we explain why Monte Carlo techniques 
cannot be used for these theories and review the Langevin equation as a numerical procedure 
for quantum field theories. In Sec. |ITT| we review the complex Langevin method and use the 
methods of supersymmetric quantum mechanics to derive the conditions that guarantee 
convergence for the expectation values. First and second-order algorithms for implementing 



the Langevin method are presented in Sec. and the results of numerical simulations are 
given and shown to be in excellent agreement with known results. Sec. contains concluding 
remarks concerning the implications of this study in regard to probability and completeness 
and proposes higher- dimensional PT-symmetric theories to which this numerical method 
could be applied. 
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II. THE LANGEVIN METHOD 



Since many PT-symmetric theories possess a real positive spectrum and a nonvanish- 
ing value for Gi, it has been speculated that PT-symmetric theories could be used to 
describe a Higgs Boson. A —gip'^/A theory is especially interesting as a theory for the 
Higgs because it has a dimensionless coupling constant and is asymptotically free |]T7 . 



However, the Schwinger- Dyson equations mentioned in Sec. |I| are too difficult to solve in 
four- dimensional space-time, where physical quantities must be calculated. Hence, a reli- 
able numerical method is needed to compute expectation values in higher dimensions. In 
this section we argue that Monte Carlo methods are ill suited for these theories, while the 
complex Langevin equation provides a robust numerical technique. 

Most often, Monte Carlo methods are used for numerical evaluations of expectation val- 
ues like those in Eq. (|1.5|) where the contour of integration is along the real axis and S{(f)) is 
real. This is achieved by choosing paths weighted according to the probability distribution, 
exp(— S'(0)). For the Hamiltonians considered in this paper, 5'(0) is complex, and neither 
a real representation nor a consistent probabilistic interpretation is known. One approach 
to use Monte Carlo is to separate the Hamiltonian into real and imaginary parts and con- 
sider y4(0) exp(— 'i/m[S'(0)]) the operator and exp(— i?e[S'(0)]) the probability distribution. 
There are two problems with this approach. Firstly, the size of /m[S'(0)] increases with the 
size of the lattice, making the numerator and denominator of Eq. ( |1.5|) very small. Con- 
sequently, numerical simulations become very difficult. Secondly, and more importantly, 
in many cases /m[S'(0)] contains more information about which paths are important than 
Re[S{(f))]. Therefore, the algorithm outlined above never samples the important paths and 
thus fails to converge to the correct answer. Specific examples of this phenomenon are given 
in Refs. |lT9|-pT|. Our studies of zero- dimensional versions of Eq. ( |1.1| ) suggest that this is 
the case for non-Hermitian but PT-symmetric theories, and that Monte Carlo methods are 
inapplicable. 

Another method for numerical calculations in quantum field theory involves the Langevin 
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equation p2|,p3 



where r is an unphysical, Langevin time, dS{(j))/d(f) gives the equations of motion for the 
Hamiltonian, and ri{T) is a stochastic variable. The function //(r) is chosen to be a real, 
Gaussian random function that satisfies the conditions 

(^(r)) = 0, {r^{TMr')) = 26{r-T'), (2.2) 

where the averaging is performed with respect to the appropriately normalized Gaussian 
probability distribution. Further, it is well known that when S{(f)) is real, the probability 
distribution, P{4>, r), associated with Eq. (pTTp is given by the Fokker-Planck equation p2|j2! 



The space-time dependence of all variables and the r dependence of 5* and (p are left imphcit 
in the above equations. These dependencies are only made explicit when relevant to a 
calculation. 

It is easy to evolve Eq. ( |2.1| ) numerically in Langevin time, r, and find expectation values. 
For real variables these expectation values are expressible as 

{0\A\0)p JD<pAi<P)P{<P,T) 



and one can show that 



(0|0)p /D0P(0,r) 



P(0,r) ^ e-^^"^), T^oo. (2.5) 



(2.4) 



Hence, 

/ D(j) A{(f))P{(j), t) J D(f) A((^)e-^('^) 



r ^ oo. (2.6) 



J D(l) P{(I),t) /D0e-^W 

That is, when S{(f)) and are both real, the physical expectation values are recovered by 
taking the unphysical, Langevin time to infinity 



More generally, an analytically continued version of Eq. ( p.5|) often holds as long as the 
supersymmetric Fokker-Planck Hamiltonian, Hpp, formed by taking dS/d(f) as the super- 
potential, has a spectrum with positive real part and a ground state that is nondegener- 
ate |T8|J2^ . When 5(0) is complex, Hpp can still have a nondegenerate ground state and 
a spectrum that is real and positive. As explained in Sec. these criteria are the correct 
ones to test for convergence of expectation values for the Hamiltonians studied in this paper. 
This is also true for several other cases studied in Refs. |25H30|. This method was successful 



in several cases including statistical mechanics problems with complex chemical potentials, 
field-theoretic calculations in Minkowski space, simulations dealing with the many Fermion 



problem [^], and even non-Hermitian Hamiltonians with complex eigenvalues |]T9 . 

The previous studies involving the complex Langevin equation focused on cases where 
the physical Hamiltonians were either Hermitian, non-Hermitian with a positive real part, 
or non-Hermitian with a real part that was negative and an imaginary part that was small. 
These cases were studied because the associated Fokker-Planck Hamiltonian, Hpp, had 
eigenvalues with a positive, real part. In contrast, PT-symmetric Hamiltonians such as 
Eq. ( |1.1| ) often have a real part that is not strictly positive and an imaginary part that 
cannot be considered small. However, in the cases we have studied, Hpp still possesses a 
spectrum that is purely real and positive. 

Reference demonstrates that the renormalized mass squared for the anharmonic 
oscillator. Ml = {Ei —Eq^, is proportional to the first nonzero eigenvalue of the associated 
Fokker-Planck Hamiltonian. This suggests that the positive real part of the spectrum of Hpp 
is a consequence of the purely real spectrum of these PT-symmetric Hamiltonians. For the 
zero-dimensional PT-symmetric Hamiltonians studied in this paper, the eigenvalues of the 
associated Hpp are always positive and real. We believe that a "PT-symmetric Hamiltonian 
with real, positive eigenvalues will always lead to a real positive part for the eigenvalues of 
Hpp. To prove this would be tantamount to finding the exact conditions necessary for a 
given PT-symmetric Hamiltonian to possess a real spectrum, and that is an open problem. 
Moreover, in Ref. |^ and elsewhere, Hermitian Hamiltonians in Minkowski space often lead 
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to an associated Hpp that is non-Hermitian. In contrast, the PT-symmetric, non-Hermitian 
Hamiltonians in this paper always lead to Fokker-Planck Hamiltonians that maintain VT 
symmetry; this is explained in Sec. |T|. 



III. COMPLEX LANGEVIN EQUATION 

To gain a deeper understanding of when the Langevin equation works and of its connec- 
tion to the Hamiltonian being studied, we begin with the complex Fokker-Planck equation. 
Allowing S{(f)) and (f) to be complex, Eq . (p . 1|) can be divided into its real and imaginary 
parts and written as two coupled equations. If one assumes the noise is purely real and uses 



the standard methods of stochastic calculus to derive Ito's formula for two variables p3 
one is lead to the complex Fokker-Planck equation. 



dr 



dS 



+ — Im 



dS 



+ TnT]P{<t>RAi:r) 



= OFp{(Pn,h)P{(^RAi;r), (3.1) 



where (pu and are the real and imaginary parts of (p respectively and S = S{(j)fi + i(f)i). 
Eq. ( |3.1|) defines a purely real probability in the complex-^ plane, but apart from a few simple 
cases [^,^, explicit constructions of P(0R,0/;r) are unknown. 



Now that we are evolving the Langevin equation in the complex plane, Eq.(2^) must be 
modified. The average over the Langevin probability must be taken as an area integral in 
the complex plane given by 

(0|^|0)p ^ / D(pR Dh A{<Pr + i<Pi)P{<Pr, <\>u r) 
(0|0)p \D4>nD4>jP{4>nAi\r) ■ 

Notice that A{(j)fi + i(f)i) is an analytic function, but P is not in general. Understanding how 
Eq.( |2.5|) and Eq.( ^.6| ) are satisfied is now much more complicated because in the limit t ^ oo 
one must show how an area integral becomes a path integral and that a real, nonanalytic 
function, P, generates the complex, analytic function exp(— 5(0)). This can be achieved in 
a formal manner by following the approach introduced in Refs. 
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For the case of ig(j)^/3 in zero dimensions with m = 0, the path integral converges when 
exp(— 5(0)) is exponentially damped. Expressing the complex variable in polar coordinates, 
(j) = rexp{i9), the Stokes' regions that are traditionally chosen for PT-symmetric theories 
are — vr < 9 < — 27r/3 and — 7r/3 < < 0, as discussed in Sec. |. These wedges are depicted 
in Fig. 1. Moreover, analytical calculations for these theories are most easily done along 
the contour where there is pure exponential damping defined by = — Svr/G and 6 = —tt/Q. 
This contour is the dashed line in Fig. 1. However, any contour whose endpoints lie in 
the appropriate Stokes wedges is acceptable. For purposes of proving convergence for the 
Langevin expectation values, the easiest contour to use is (p = (pR — ib, where b is any 
finite constant. Along this contour, exp{—S{(f))) = exp{—ig(j)^/3) ~ {oscillatory term) * 
exp(— 5f0|,6), and is therefore damped as (pR ±oo. This contour is the solid line in Fig. 1. 

The large r behavior of the expectation values given by Eq. ( p.2|) is discovered by shifting 
integration variables, 0/ — > 0/ — b. This shift does not affect the endpoints of integration 
and has a Jacobian of one. Any analytic function A{(j)ii — ib + icpi) can be Taylor expanded 
about the contour, (pc = (pR ~ i^- This allows us to express the expectation values as 
(0|A|0)p _ ID<PrD<Pi [e'^A{<Pc))P{<PR,<Pi-b;r) 



(0|0)p JD(PrD<Pj Pi(PR,(Pi-b;T) 

where 



(3.3) 



X = h4r- (3-4) 



Integrating Eq. (|3.3| ) by parts infinitely many times gives: 

(0|A|0)p _ !D<Pr A{<Pc)Peff{<Pc.r) 
(0|0)p JD<PRP,jf{<Pc,r) ^ 

where 



(3.5) 



P. 



J D<Pje-'^P{<PR,<Pj-b;r). (3.6) 



We assume that P vanishes at infinity rapidly enough so that all of the boundary terms 
from the integration by parts are zero. (In the denominator of Eq. (|3.3| ) exp(— zx) can be 
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introduced for free because all but the zeroth order term in x integrate to zero.) Note that 
P(.ff is an analytic function of <pc = 'Pr ~ ib, not a function of (pR and b separately. We 
see this by using exp{—ix)P{(t>R, (pi — b;T) = Pi^cpR — icpi, (pi — b; r) and then shifting the 
integration variable (pi (pi + b, so that Pg// = / D(pi P{(pc — i(pi,(pii'T)- As a result, 
Eq.(|3.5D can be equivalently written as 

(0|^|0)P ^ I-o^%bD<Pc A{(Pc)Peff{<Pc,T) 

(0|0)P r--%D(Pc Peffi<Pc,r) • 

We now derive a pseudo Fokker-Planck equation for Peff{(pc,T). From Eqs. (|3.6| ) and 



dPeffi^c, r) _ r ^ .i^dP{4>R, (pi - b; 



D(pi e 



T 



dr J dr 

J D(Pi 0'Jle~'^P{(PR, (Pi -b-r), (3.8) 



where O'pp = exp{—ix)OFp{(pR, (pi — b) exp(ix). Using the relations 



d(pR d(pR 
d ■ d d 



0({)i 0(pi 0(pR 

e-'^F{(Pc + icpM"" = F{<Pc), (3.9) 
it is straightforward to show that 

O'Jl = 3^ f + + + (3.10) 

d(pR \d(pR d<pR / d<pi d(p 



The last term of Eq. ( p.lO|) is a total derivative in 0/, and therefore, disappears from the 



right side of Eq. (p.8|) , again assuming P vanishes rapidly at infinity. Since the remaining 
terms of O^Jp do not depend on 0/, they can be pulled out in front of the integral over 
D(pi. Using the fact that d/d(pR = d/d(pc on an analytic function of (pc, Eq. (pl8|) becomes 
Eq. (pD with <P = <Pc and P(0, r) = Pe//(0c, r) 



dr 
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That is, there is a pseudo-Fokker-Planck equation that defines a complex analytic function, 
Peff, that is just the analytic continuation of the Fokker-Planck equation for real variables. 

For cases where > 3 in Eg. ( |1 . 1| ) , a similar derivation gives the same result. The only 
subtlety is in choosing the correct contour. For any a contour with finite endpoints, 
(—6, ~b cot 9l) and (—6, 6 cot 6*/?), as defined by Eq.( |T]^), can be deformed into the contour 
0c = 0i? ~ "ib as shown in Fig. 2. Consequently, using the methods above, an area integral 
over the strip —bcotOi < (pR < b cot Or, — oo < 0/ < oo is expressible as a path integral 
over the contour 0c = 0_r — ib plus boundary terms involving derivatives of P. As b grows 
larger, the area integral approaches an integral over the entire complex plane. We expect 
the boundary terms to approach zero because the probability of finding a particle at ±oo 
should go to zero. This is seen in Fig. 3 (Sec. Following the derivation for ig(j)^/3, we 
are again led to Eq.( p.ll| ) 



Thus, the problem of understanding the r — > cx) behavior of expectation values has been 
reduced to one that is formally identical to that for real variables. We now use the methods of 
Parisi and Sourlas |]53| who first discovered the hidden supersymmetry in classical stochastic 
equations. 

If we express Eq.( ^.ll| ) in terms of p{(pc,T) = -Pe//(0C) t) exp(S'(0)/2), we obtain the 
Schrodinger equation 

-8p(4'c,r) „ ,_( S ldS\(d ias\ 

As claimed, Hpp is the supersymmetric Hamiltonian formed from the superpotential 
dSjdcpc- Since S is VT symmetric, dS/dcpc is anti-PT symmetric, and Hpp is VT- 
symmetric, as claimed in the last section |^|. Expanding Eq. (|3.12|) yields 



If the time-independent version of Eq. ( |3.12| ), 



fffpH-rWc) = At*r('#'c), (3.14) 
is well posed and the eigenfunctions "^^^{(pc) are complete, then 
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P(0c, r) = E afeV^r(0c)e-''^-^ (3.15) 

fc=0 

Note that \1/q'^(0(7) = exp(— S'(0c')/2) is an eigenfunction of Hpp with Aq = 0. Therefore, 
Eq. ( |3.15|) becomes 



p(0c, r) = Ce-(^/^)^(*-) + E «.^r(0c)e-^^-. (3.16) 

k=l 

Moreover, if the spectrum of Hpp is such that 

Re[\k] > 0, A; > 0, (3.17) 

it follows that 

r ^ oo. (3.18) 

The T dependence of p{(f)c,T) has disappeared in this limit. That implies dPeff/dr = 
and signals that the system has reached equilibrium. Expressing Pf.ff{(j)ciT) in terms of 
p{4>Cj 't) and taking the limit r ^ oo gives analytically continued versions of Eq. ( |2.5| ), and 
thus Eq. (|2.6| ), in terms of Peff- As a result, Langevin expectation values are shown to 
converge to the right side of Eq. ( |1.5| ) as r ^ oo. This result is true for our zero-dimensional 
PT-symmetric theories as long as the ground state is nondegenerate. There is no evidence 
that PT-symmetric theories possess a degenerate ground state, so for the purposes of this 
paper, we will not consider this a possibility. 

Thus, if the supersymmetric Hamiltonian Hpp formed from the superpotential, dS/d(f)Ci 
has a spectrum satisfying Eq. (|3.17|) , the Langevin method should work as a calculational 



procedure. Explicitly, we have shown that analytic continuations of Eq. (p.5|) and Eq. ( p.6|) 
with P replaced by Peffi 

Peffi^c, r) ^ e-^(<^^), r ^ oo (3.19) 



and 



/ D^R D^j A{<Pp + t<Pj)P{<PR, h\ r) _ r-%b D<Pc A{^c)Peff{<Pc, r 



J D(t>R D<l>j P(0^, 0,; r) /-J^^, D^c Peffi^c, r) 



r ^ oo (3.20) 
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follow if Hpp has a nondegenerate ground state, wave functions that are complete, and a 
spectrum with positive real part. 



IV. NUMERICAL METHODS AND RESULTS 

In this section we apply the complex Langevin method to massless versions of Eq. ( |1 . 1| ) in 
zero and one dimension and calculate the same time one-point and two-point disconnected 
Green's functions for the cases = 3 and = 4. We begin by proving (under certain 
assumptions) that Eq. (|3.17|) holds in zero dimensions and explaining the algorithms we 



have used to implement simulations. We then extend these results to their one-dimensional 
analogues. 

A. Zero-Dimensional Theories 



A recent paper by Dorey et al. |16| shows that PT-symmetric Hamiltonians of the form, 



n = -(90)" - (z0)"'" - a{i(l)f'-\ (4.1) 

where M and a are real and boundary conditions have been chosen as in Sec. |, have a real 
positive spectrum if the conditions a < M and M > 1 are both satisfied (We have set 
/ = in the Hamiltonian given by Dorey et al..) This proves (for a = 0, M = N/2) that 
massless versions of Eq. ( |1 . 1|) have a positive real spectrum. 

In zero-dimensional studies of Eq. ( |1 . 1|) with m = 0, Eq. ( p.l3|) gives 



- 2 ~ ^ 

where the contour, (pc, is within the Stokes' wedges explained in Sec. | and used by Dorey 
et al. Making the change of variables (p {2/g)^^^(f), Eq. ( |4.2| ) becomes Eq. ( [4.1| ) with 
a = — 1 and M = — 1. Thus, a = M, and as explained in Ref. [0, this implies that 
the Hpp given in Eq. ([4.2|) has one zero eigenvalue, which we have already demonstrated. 
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and that all of the remaining eigenvalues are real and positive. Consequently, Eq. (p.l7|) 
holds for N >2, and this implies that the complex Langevin method will work. 

Further, the large behavior of the wave functions for the eigenvalue problems defined 
by Eq. ( ^.2| ) must have the form 




by WKB. Apart from a factor of two, the asymptotic form of the wave functions have 
exactly the same form as exp(— S'(0)) = exp{g{i(f))^ /N). Thus, the Stokes' wedges that de- 
fine regions of convergence for the path integral are exactly the same as those that demand 
\l/(±oo) = 0. This is equivalent to noting Eq. ( |1.4|) with — > 2{N — 1) is identical to 
Eq. (|1.2|) . That is, the wedges of convergence for the path integrals defined by exp(— S'((/))) 
are preserved by the boundary conditions for the wave functions of the Fokker-Planck Hamil- 
tonian. 

The most general form of the Langevin equation is 

|^ = F(0(r))+r/(r) (4.4) 

The simplest discretization of this is Euler's method (a first-order algorithm) and is explicitly 
given by 

0(j + 1) = 0(j) + /i(F(0(j)) + vU)) = m + e'Fi<Pij)) + eV(j), (4.5) 

where j is an index for a Langevin-time step, h is the spacing in Langevin time, e 
and the h dependence of ri{T) has been made explicit 

V'U) = evij), iv'ijWik)) = 2S,,,. (4.6) 

This form for rj^r) follows from the the normalization of the Gaussian probability distri- 
bution because Eq. (|2.2| ) implies ?7^(r) ~ 5(0) ~ on a lattice. For this and the second 
order algorithm given below, there are numerical instabilities for large values of 0. The worst 
instabilities arise when the potential is —g(j)'^/4. For —g(f)'^/4, it was necessary to restrict 
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the absolute magnitude of in order to avoid these instabihties. A typical plot of the path 
followed by (f) in the complex plane is shown in Fig. 3. 

One must then take e ^ 0. Limiting values were obtained by fitting the data with second 
degree and third degree polynomials in e. These fits are similar to those seen in Fig. ^, which 
is for the one- dimensional case. Errors are calculated by collecting the simulation data in 
bins of a given size and computing the standard deviation of the means of the bins. The 
maximum error as a function of bin size is taken to be the error for the simulation. In Table | 
the numerical results obtained using Euler's method are compared with exact values given 



in Ref. ^^j- The parts of the one-point and two-point disconnected Green's functions that 
are known to vanish (e.g. Re[Gi\) have errors larger than their magnitude as e — 0. 

Euler's method is expected to converge linearly in e as e ^ 0. Therefore, a more accurate 
second-order in e method is desirable. A second-order Runge-Kutta algorithm that led to 



good results in previous studies and was first developed in Ref. |35 



IS 



0(j + 1) = 0(j) + \e\F{<P{3)) + + eV(j). (4.7) 

In our studies this method is more stable numerically than Euler's method, and therefore, 
allows the inclusion of more data. Limiting values were obtained by fitting the data with 
second degree and third degree polynomials in e with the linear term set equal to zero. 
(These fits are similar to those seen in Fig. ^ for the one-dimensional case.) The results 
obtained using this algorithm are compared with exact values and the result of Euler's 
method in Table |I[ There is good agreement. 

It should be noted that for the case = 4, 0(0) has to be chosen in the lower half of 
the complex-0 plane or else the numerical simulations are unstable. This is in accord with 
the WKB wedges needed to properly define the boundary conditions as explained in Sec. |. 
This restriction also holds in one dimension for the initial configuration at r = 0. 
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B. One-Dimensional Theories 



In this section we apply the complex Langevin method to massless versions of Eq. { \L . 1|) 
in one dimension with = 3 and = 4. The potentials are the same as those discussed 
in the last section, but now there is a second derivative in physical time, —d'^(f){T,t)/dt'^, 
present in all of the equations, and //(r) becomes t dependent. Thus, the eigenvalue problem 
for Hpp becomes a partial differential equation. The spectra of eigenvalue problems for VT- 
symmetric partial differential equations has not been studied, and we do not have a proof 
that they are real. But, following the examples of zero- dimensional theories, we assume 
the spectrum of Hpp has a positive real part. The results of our simulations support this 
assumption. 

These one-dimensional theories are more computationally expensive because they give 
rise to Langevin equations that are partial differential equations. Consequently, the second- 
order algorithm in the last section is useful. For numerical simulations of one-dimensional 
theories, we have to introduce a physical time lattice and — 9^0(r, t)/(9t^ becomes 

_ (0(j,/ + l)-20(j,O + 0(j,/-l)) ^ 
a? 

where / is an index for real time and a is the spacing in physical time. The algorithms used 
in the previous section are applicable here as well, but their form has changed slightly. The 
algorithms now contain Eq. ( [4.8|) as part of F and the explicit lattice dependence of f]{T,t) 
is such that 

r]'{j,l) = Vh^vUJ), {v'{j,l)r]'{k,m)) = 25,- ^5^,^. (4.9) 

Further, in these algorithms there is always an associated with Eq. ( [4.8|) , and thus, the 
simulation is unstable unless e << a. However, there are no instabilities similar to those 
encountered for —gip'^/A in the zero dimensional case. For fixed values of a, we compute at 
various values of e and take the limit e ^ 0, giving the expectation values as a function 
of a. A typical fit for this process is shown in Fig. 4. We then take a — ^ and obtain 
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the expectation values in the continuum hmit. A fit used to extrapolate the value of —iGi 
for the potential —g(f)'^/A is shown in Fig. 5. In Table y the numerical results for the 
expectation values are compared with numerically integrated quantum-mechanical values 
given in Ref . . We have chosen g such that Eq. ( |1 . 1| ) corresponds with the Hamiltonians 



in Ref. |M. The agreement of our results with previous work is excellent. 



V. CONCLUSIONS AND SPECULATIONS 



This paper and Ref. have shown that PT-symmetric theories are amenable to the 
methods of quantum field theory. The previously-used Schwinger- Dyson method is very ac- 
curate but is difficult to apply to higher-dimensional theories. Here we show how a numerical 
method based on the complex Langevin equation can be used to obtain precise results. We 
believe that this numerical method can be applied to higher-dimensional theories and plan 
to use it for future calculations. This work represents an important step towards a test of 
the physical applicability of PT-symmetric Hamiltonians. 

An interesting implication of this study concerns the probability and completeness of 
PT-symmetric theories. The argument for the success of the complex Langevin method 
in Sec. |Tl| crucially depends upon the eigenfunctions of Hpp being complete ( p.l5| ). Simi- 
larly, extremely accurate results in Ref. [0 crucially depend on the completeness of VT- 
symmetric eigenfunctions. The possible connection between the eigenvalues of the Fokker- 
Planck Hamiltonians and the Hamiltonians being simulated suggests a possible connection 
between the eigenfunctions as well. Perhaps proving completeness for one of these sets of 
eigenfunctions would imply the completeness of the other set. 

Perhaps the most important implication for PT-symmetric theories is that there is an 
implicit probability distribution defined by Eq. and seen in Fig. 3. Moreover, this study 
suggests that expectation values must be interpreted as area integrals of observables weighted 
by this real probability distribution. Previous studies of the zeros of eigenfunctions of VT- 
symmetric theories suggest that the zeros interlace and become dense in a narrow band in the 
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complex-^ plane |3g]. Hence, it was conjectured that completeness may have to be defined 
in terms of area integrals. We suspect a strong connection between the two results and 
speculate that a consistent probabilistic formulation of PT-symmetric theories can only be 
achieved in terms of area integrals. A probabilistic interpretation of "PT-symmetric theories 
would be a major advance. 

A recent study of the complex Langevin equation demonstrates that problems with 
convergence arise when expectation values are complex but the fixed points of the Langevin 
equation lie on the real axis, and as a result, the field spends most of its time on the real 



axis and away from the correct average value It was demonstrated there that this 

problem could be avoided by moving the fixed points into the complex plane. In the present 
study the expectation values are pure real or pure imaginary. The fixed point of the Langevin 
equation is zero, which is on the real axis, but the simulations converged without moving the 
fixed point into the complex plane. This seems to be because the path of the deterministic 
equation (r/ = 0) is given by 

0(r) = (5.1) 

[{N -2){t + C)]—^ 

where C is an arbitrary constant given by the initial condition. The field is attracted to the 
imaginary axis even though the fixed point is on the real axis. Thus, it appears to be crucial 
that the deterministic path to or between the fixed points is in the complex plane, but not 
that the fixed points themselves lie in the complex plane. The simulations for PT-symmetric 
theories work in a straightforward manner; the fixed points do not need to be adjusted. We 
believe that the complex Langevin equation works so nicely for PT-symmetric theories 
because these theories are inherently complex. As a result, it seems that PT-symmetric 
theories may provide a class of toy models with new and interesting properties that are 
especially well suited for probing the intricacies of the complex Langevin equation. 

The most exciting proposed application of these non-Hermitian theories is to Higgs the- 
ories. We believe that the numerical method presented in this paper should allow one to 
compute the mass of the Higgs particle in four- dimensional versions of these theories. Before 
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this can be done, however, renormahzation must be thoroughly understood. Studies of field 
theories in fewer than four physical dimensions represent a simpler alternative. In particular, 
exact results for the scaling exponents of an igcj)^ /?> theory in two and three dimensions can 
be obtained by relating them to the Lee- Yang edge singularity in one and two dimensions 
respectively [^,^. Evaluation of these exponents by the methods introduced here are in 
progress. 



20 



TABLES 



TABLE I. Numerically determined values of id = i(0|(^|0)/(0|0) and -G2 = -(0|(/)2|0)/(0|0) 
using Euler's method and a second-order method for zero-dimensional igcj)^ /3 and —gcf)'^ /A theories 
with g = 1/2. These limiting values were determined by fitting the simulated data to polynomials 
of second and third degree in e. For the second-order algorithm, the term linear in e is set to zero. 
Exact results are listed to four significant digits in the first column. Note that the values listed for 
the second-order algorithm are indeed more precise than the results using Euler's method. 





N iGf^^ 




^^2ndorder 


^exact 


_^Euler 


^2ndorder 


3 0.9185 

4 1.163 


0.9198(14) 

1.166(3) 


0.9194(7) 

1.164(1) 




0.9560 


0.9623(51) 


0.9595(14) 



21 



TABLE II. Numerically determined values of iGi = i(0|(/>|0)/(0|0) and -G2 = -(0|(/>2|0)/(0|0) 
using Euler's method and a second-order method for one-dimensional —g{i(j))^/N theories with 

= 3 and = 4. In order to compare our results with Lagrangians used in previous studies, we 
have set g = N/2. These continuum limit values were determined by fitting the simulated data to 
third degree polynomials in e for a given a (for the second-order algorithm, the term linear in e is 
set to zero), and then fitting the results of those fits with second degree polynomials in a?. Exact 
results obtained by direct numerical integration of the quantum mechanical problem are listed to 
four significant digits in the first column. The exact value of —G2 for A^ = 4 has never been 



directly calculated, so instead we list the results of the variational calculation given in Ref. |14|. 
Note that the values listed for the second-order algorithm are indeed more precise than the results 
using Euler's method. 





N 


^^cxact 


^^Eulcr 


^^2ndordcr 


^var 


_ ^Eulcr 


^2ndordcr 


3 
4 


0.5901 
0.8669 


0.5890(5) 
0.8654(6) 


0.5898(2) 
0.8670(3) 




0.5182 


0.5171(9) 


0.5183(4) 
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FIGURES 




FIG. 1. Stokes' wedges for a massless, zero-dimensional igcf)^ /2 theory. The shaded areas are the 
regions within which the path integral for this theory converges because exp(— S'(0)) is exponentially 
damped. The usual path of integration is represented by the dashed line and extends from — oo 
to the origin along the ray defined by Or = — 57r/6 and then from the origin to oo along the ray 
defined by 0l = — 7r/6. Along this contour exp(— S'(^)) is purely exponentially damped. In order to 
prove that the Langevin expectation values have the desired behavior as r ^ oo, it is easiest to use 
the smooth contour represented by the solid line. The solid line contour is defined as ^ = — ib 
and extends from —oo — ib ^ oo — ib. 
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FIG. 2. Typical example of Stokes' wedges for a massless, zero-dimensional —g{i(t))^/N theory. 
The shaded areas are the regions within which the path integral for a given N converges because 
exp(— S'((/>)) is exponentially damped. The usual path of integration, even for a finite contour, is 
represented by the dashed line and extends from (— b,—b cot 9 l) to the origin and then from the 



origin to (—6, 6cot 0^), where 6l and 0r are given in Eq. 1.2. In order to prove that the Langevin 
expectation values have the desired behavior as r ^ cxd, it is most convenient to use the smooth 
contour represented by the solid line. The solid line contour is defined as (p = (pR — ib and extends 
from (— 6, — ^cot^i) (— 6, 6cot 
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FIG. 3. Plot of the field in the complex plane for the potential —g^i^/A with e = 0.3 using 
the second-order algorithm. Each point corresponds to the value of (j) for a value of the fictitious 
time r. The absolute magnitude of <j) was restricted to be less than 19 in order avoid numerical 
instabilities. (For simulations of —g<^'^/4 in one dimension, this restriction was unnecessary.) The 
first 10,000 points are plotted. The field started at the point (0.5, —0.1) and then followed a path 
towards the imaginary axis. It then traveled from side to side forming a cloud of points that 
averaged to the value —1.1687. 
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FIG. 4. Plot of values for —iGi for a real time spacing of a = 0.5 versus the square of the 
fictitious time spacing, e^, for a —g4)^/A potential in one dimension. The solid line is a fit that is a 
linear plus quadratic in e for values computed using Euler's method (diamonds). The dashed line 
is a cubic fit in e, without a term that is linear in e, for values computed using the second-order 
algorithm (squares). Similar fits were performed for each value of a shown in Fig. 5. 
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FIG. 5. Extrapolated values of —iGi versus real time spacing for a —g<p^/A potential in one 
dimension. Lattice expectation values for these theories only depend upon even powers of a, so 
fits are done in powers of o^. Both fits shown here are linear plus quadratic in a^. Diamonds show 
values computed using Euler's method, and the solid line is a fit for them. Squares show the results 
from the second-order algorithm; the dashed line is the fit to them. The continuum limit values 
obtained from these plots are given in Tab. 2. Similar fits were performed for all of the data in 
Table 2. 
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